/*  README:
*       Algoritmo de Gauss Legendre(Sequencial)
*
*           Algoritmo que converge na ordem de 2 casas decimais com precisao de 10.000.000 casas decimais,
*           portanto temos 25 iteracoes para seu calculo na main;
*           para alterar a quantidade de iteracoes deve-se alterar o iterador de 25 para o desejado na funcao;
*           o tempo eh fornecido em milisegundos; defina a precisao na variavel "prec"
*
*/

#include<stdio.h>
#include<stdlib.h>
#include<gmp.h>
#include<time.h>

//Inicializa a0 com 1
void inicializaA0(mpf_t *a0){
    mpf_set_str(*a0, "1",10);
}

//Inicializa b0 com 1/sqrt(2)
void inicializaB0(mpf_t *b0){
    mpf_t aux;
    mpf_init(aux);

    mpf_set_str(*b0, "1",10);
    mpf_set_str(aux, "2",10);

    mpf_sqrt(aux,aux);

    mpf_div(*b0, *b0, aux);

    mpf_clear(aux);
}

//Inicializa t0 com 1/4
void inicializaT0(mpf_t *t0){
    mpf_t aux;
    mpf_init(aux);

    mpf_set_str(*t0, "1",10);
    mpf_set_str(aux, "4",10);

    mpf_div(*t0, *t0, aux);

    mpf_clear(aux);
}

//Inicializa p0 com 1
void inicializaP0(mpf_t *p0){
    mpf_set_str(*p0, "1",10);
}

void proximoA0(mpf_t *a, mpf_t b){
    mpf_t aux;
    mpf_t aux2;
    mpf_init(aux);
    mpf_init(aux2);

    mpf_add(aux2, *a, b);
    mpf_set_str(aux, "2",10);

    mpf_div(*a, aux2, aux);

    mpf_clear(aux);
    mpf_clear(aux2);
}

void proximoB0(mpf_t *b, mpf_t a){
    mpf_t aux;
    mpf_init(aux);

    mpf_mul(aux, a, *b);
    mpf_sqrt(*b, aux);

    mpf_clear(aux);
}

void proximoT0(mpf_t *t, mpf_t p, mpf_t a0, mpf_t a1){
    mpf_t aux;
    mpf_init(aux);

    mpf_sub(aux, a1, a0);
    mpf_pow_ui(aux, aux, 2);
    mpf_mul(aux, aux, p);

    mpf_sub(*t, *t, aux);

    mpf_clear(aux);
}

void proximoP0(mpf_t *p){
    mpf_t aux;
    mpf_init(aux);

    mpf_set_str(aux, "2", 10);

    mpf_mul(*p, *p, aux);
}

void calculaPi(mpf_t *pi, mpf_t a, mpf_t b, mpf_t t){
    mpf_add(a, a, b);
    mpf_pow_ui(*pi, a, 2);

    mpf_set_str(b, "4", 10);
    mpf_mul(t, t, b);

    mpf_div(*pi, *pi, t);
}

int main(){
    clock_t tInicio, tFim;
    tInicio = clock();

    int i, prec = 33211928;//10.000.000 * log(10)/log(2)
    mpf_t a0;
    mpf_t b0;
    mpf_t t0;
    mpf_t p0;
    mpf_t aAux;
    mpf_t bAux;
    mpf_t tAux;
    mpf_t pi;
    mpf_t piAnterior;

    mpf_set_default_prec(prec);

    mpf_init(a0);
    mpf_init(b0);
    mpf_init(t0);
    mpf_init(p0);
    mpf_init(aAux);
    mpf_init(bAux);
    mpf_init(tAux);
    mpf_init(pi);
    mpf_init(piAnterior);

    inicializaA0(&a0);
    inicializaB0(&b0);
    inicializaT0(&t0);
    inicializaP0(&p0);

    for(i=0; i<25; i++){
        mpf_init_set(aAux, a0);
        mpf_init_set(bAux, b0);
        mpf_init_set(tAux, t0);

        proximoA0(&a0, b0);
        proximoB0(&b0, aAux);
        proximoT0(&t0, p0, aAux, a0);
        proximoP0(&p0);

        calculaPi(&pi, aAux, bAux, tAux);

        gmp_printf("Interacao %d: PI = %.30Ff \n",i+1,pi);
    }

    tFim = clock();
    printf("\nNumero Total de Interacoes: %d    Tempo Decorrido:%ld ms\n\n", i,(tFim-tInicio)/(CLOCKS_PER_SEC/1000));

    return 0;
}
